---
title: "AL Reliability Analysis_2023-07-17"
author: "Sammy"
date: "7/17/2023"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load libraries, include=FALSE}
library(tidyverse) 
library(irr) # interrater reliability
library(irrCAC) # gwet's AC
```

```{r load data, include=FALSE}
## dataframes ##
# dataframes are set up so each row is possible AL that was rated/denoted by at least on rater. There is a column that donates the specific tooth/the tooth code, what area the AL was marked in (prenatal enamel, postnatal enamel, or dentine), and a column for each individual or team of raters that contains their rating (RL, MAL, or HAL). 

## dentine_rel --> dentine reliability df
## mal_hal_rel --> MAL and HAL reliability df 
## overall_rel --> Presence and absence of AL reliability df
## prenat_rel_full --> Prenatal reliability df
## postnat_rel_full --> Postnatal reliability df
## rating_cols --> vector of colnames/rater names
## rel_data_raw --> original df before cleaning 

## vectors ## 
## jr --> junior-level raters
## senior --> senior-level raters
## mid --> mid-level raters
## no_al --> raters with no AL experience
## yes_al --> raters with AL experience 
## no_decid --> raters with no deciduous teeth experience 
## yes_decid --> raters with deciduous teeth experience 
## under_200 --> raters with less than 200 slides
## up_500 --> raters with 200 - 500 slides
## plus_500 --> raters with 500+ slides 
## rating_cols --> column names/names of raters 
```

```{r lights weighted function, include=FALSE}
lights_weighted <- function(raters, weighting, ratings){
weights_table <- data.frame(grouping = character(0), kappa = numeric(0)) 
for (i in raters){
  for (j in raters){
    if (i == j){} # do not run rater against self
    else if (paste0(j,"+",i) %in% weights_table$grouping){} # do not repeat pairing 
    else{
   x <- kappa2(ratings[c(i,j)], weight = weighting)
   value <- x$value
   temp <- data.frame(grouping = paste0(i,"+",j), kappa = value)
   weights_table <- rbind(weights_table, temp)}
  }
}
light_kappa = round(mean(weights_table$kappa),3)
return(light_kappa)
}


## make ratings ordered factors 
mal_hal_rel[rating_cols] <- lapply(mal_hal_rel[rating_cols], factor, levels = c("retzius", "mal", "hal"), ordered = TRUE)

prenat_rel_full[rating_cols] <- lapply(prenat_rel_full[rating_cols], factor, levels = c("retzius", "mal", "hal"), ordered = TRUE)

postnat_rel_full[rating_cols] <- lapply(postnat_rel_full[rating_cols], factor, levels = c("retzius", "mal", "hal"), ordered = TRUE)

## export function
#saveRDS(lights_weighted, file = "path")
```

```{r groupings table, include=FALSE}
groups_table_df <- data.frame(group = c("Senior", "Mid", "Junior", "No Deciduous Teeth Experience", "Deciduous Teeth Experience", "No AL Measurement Experience Group", "AL Measurement Experience Group", "< 200 Slides", "200 -500 Slides", "500+  Slides"),
                rater3 = c(NA,NA,"X","X",NA,"X",NA,"X",NA,NA),
                rater5 = c(NA,NA,"X","X",NA,"X",NA,"X",NA,NA),
                rater6 = c(NA,NA,"X","X",NA,"X",NA,NA,"X",NA),
                rater8 = c(NA,"X",NA,NA,"X",NA,"X",NA,"X",NA),
                rater4 = c(NA,"X",NA,NA,"X",NA,"X",NA,"X",NA),
                rater2 = c(NA,"X",NA,"X",NA,NA,"X","X",NA,NA),
                rater1 = c("X",NA,NA,"X",NA,"X",NA,NA,NA,"X"),
                rater7 = c("X",NA,NA,NA,"X",NA,"X",NA,NA,"X"))

library(htmlTable)
groups_table <- groups_table_df %>%
  addHtmlTableStyle(col.rgroup = c("none", "#F7F7F7"),
                    col.columns = c("none", "#F1F0FA")) %>% htmlTable(rnames = FALSE)
```

# "NNL Doubling" 
```{r notes and removing NNL doubling, include=FALSE}
# number of instances of marking NNL doubling 
doubling <- rel_data_raw %>% filter(str_detect(note_on_line_characteristics, regex("doub", ignore_case = T)))
n_doubling <- length(doubling$line_code) # 3 cases 

# number of cracks 
cracks <- rel_data_raw %>% filter(str_detect(note_on_line_characteristics, regex("crack", ignore_case = T)))
n_cracks <- length(cracks$line_code)
```

```{r prevalance table, include=FALSE}
## overall rel (yes/no AL)
length(which(overall_rel[,rating_cols] == "no"))/(length(rating_cols)*nrow(overall_rel)) 
length(which(overall_rel[,rating_cols] == "yes"))/(length(rating_cols)*nrow(overall_rel)) 

## MAL vs HAL vs Retz reliable --> removes dentine 
length(which(mal_hal_rel[,rating_cols] == "hal"))/(length(rating_cols)*nrow(mal_hal_rel)) 
length(which(mal_hal_rel[,rating_cols] == "mal"))/(length(rating_cols)*nrow(mal_hal_rel)) 
length(which(mal_hal_rel[,rating_cols] == "retzius"))/(length(rating_cols)*nrow(mal_hal_rel))  

## long format
mal_hal_rel <- as.data.frame(mal_hal_rel)
mal_hal_long <- mal_hal_rel %>% pivot_longer(cols = all_of(rating_cols),
                             names_to = "rater",
                             values_to = "rating")

## mal_hal_retz ratings table
table(mal_hal_long$rater, mal_hal_long$rating)
```

### Overall Reliability
```{r overall reliability results, include=FALSE}
# counts 
apply(overall_rel[rating_cols],2,table) 

# light's kappa due to more than 2 raters
overall_rel_result <- kappam.light(overall_rel[rating_cols])

## Gwets 
overall_rel_gwet <- gwet.ac1.raw(overall_rel[rating_cols]) 
```

### HAL vs. MAL vs. Retzius Line 
```{r mal vs. hal vs. retz reliability results, include=FALSE}
# counts 
apply(mal_hal_rel[rating_cols],2,table) 

# lights kappa due to more than 2 raters
mal_hal_rel_result <- kappam.light(mal_hal_rel[rating_cols])

# gwet's AC to check
mal_hal_rel_gwet <- gwet.ac1.raw(mal_hal_rel[rating_cols]) 

## weighted 
# lights 
mal_hal_weight <- lights_weighted(raters = rating_cols, weighting = c(0,0.5,1), ratings = mal_hal_rel)

# gwets 
mal_hal_weight_gwets <- gwet.ac1.raw(mal_hal_rel[rating_cols], weights = "linear")

mal_hal_unw_gwets <- gwet.ac1.raw(mal_hal_rel[rating_cols], weights = "unweighted")
```


## Prenatal Reliability   
```{r prenatal reliability results, include=FALSE}
# counts 
apply(prenat_rel[rating_cols],2,table)

# lights kappa due to more than 2 raters
prenat_rel_result <- kappam.light(prenat_rel[rating_cols])

# gwet's AC to check
prenat_rel_gwet <- gwet.ac1.raw(prenat_rel[rating_cols])

## weighted kappa for mal/hal ratings
# lights 
prenat_weight <- lights_weighted(raters = rating_cols, weighting = c(0,0.5,1), ratings = prenat_rel_full)

# gwets's -- use bc so many less HAL 
prenat_weight_gwets <- gwet.ac1.raw(prenat_rel_full[rating_cols], weights = "linear") 

prenat_unw_gwets <- gwet.ac1.raw(prenat_rel_full[rating_cols], weights = "unweighted") 
```

## Postnatal Reliability    
```{r postnatal reliability results, include=FALSE}
# counts 
apply(postnat_rel[rating_cols],2,table)

# lights kappa due to more than 2 raters
postnat_rel_result <- kappam.light(postnat_rel[rating_cols])

# gwet's AC to check
postnat_rel_gwet <- gwet.ac1.raw(postnat_rel[rating_cols])

### weighted kappas
# lights
postnat_weight <- lights_weighted(raters = rating_cols, weighting = c(0,0.5,1), ratings = postnat_rel_full) 

# gwets
postnat_weight_gwets <- gwet.ac1.raw(postnat_rel_full[rating_cols], weights = "linear")

postnat_unw_gwets <- gwet.ac1.raw(postnat_rel_full[rating_cols], weights = "unweighted")
```

## Dentine Reliability    
```{r dentine reliability results, include=FALSE}
# counts 
apply(dentine_rel[rating_cols],2,table)

# lights kappa due to more than 2 raters
dentine_rel_result <- kappam.light(dentine_rel[rating_cols])

# gwet's AC to check
dent_rel_gwet <- gwet.ac1.raw(dentine_rel[rating_cols])
```

```{r variation in mal hal retz categorization, include=FALSE}
## check number/percent of cases where it was just a MAL vs Retz decision, just MAL vs HAL decision, or if it was a toss up between all three ratings 
mal_hal_rel <- mal_hal_rel %>%
 rowwise() %>%
 mutate(mal_yn = +any(str_detect(c_across(rating_cols), "mal"))) %>%
 mutate(hal_yn = +any(str_detect(c_across(rating_cols), "hal"))) %>%
 mutate(retz_yn = +any(str_detect(c_across(rating_cols), "retz"))) %>%
 mutate(all_mal = ifelse((mal_yn == 1 & hal_yn == 0 & retz_yn == 0), 1, 0)) %>% 
 mutate(all_hal = ifelse((mal_yn == 0 & hal_yn == 1 & retz_yn == 0), 1, 0)) %>% 
 mutate(mal_hal = ifelse((mal_yn == 1 & hal_yn == 1 & retz_yn == 0), 1, 0)) %>% 
 mutate(mal_retz = ifelse((mal_yn == 1 & hal_yn == 0 & retz_yn == 1), 1, 0)) %>%
 mutate(hal_retz = ifelse((mal_yn == 0 & hal_yn == 1 & retz_yn == 1), 1, 0)) %>%
 mutate(mal_hal_retz = ifelse((mal_yn == 1 & hal_yn == 1 & retz_yn == 1), 1, 0))

# counts of each category per line
mal_hal_rel <- mal_hal_rel %>%
  mutate(count_mal = str_count(rater3, "mal") + str_count(rater6, "mal") + str_count(rater5, "mal") + str_count(rater4, "mal") + str_count(rater8, "mal") + str_count(rater2, "mal") + str_count(rater1, "mal") + str_count(rater7, "mal")) %>% 
  mutate(count_hal = str_count(rater3, "hal") + str_count(rater6, "hal") + str_count(rater5, "hal") + str_count(rater4, "hal") + str_count(rater8, "hal") + str_count(rater2, "hal") + str_count(rater1, "hal") + str_count(rater7, "hal")) %>%
  mutate(count_retzius = str_count(rater3, "retzius") + str_count(rater6, "retzius") + str_count(rater5, "retzius") + str_count(rater4, "retzius") + str_count(rater8, "retzius") + str_count(rater2, "retzius") + str_count(rater1, "retzius") + str_count(rater7, "retzius"))
  
# separate table for counts 
line_cat_table <- data.frame(combo = c("all_mal","all_hal","mal_retz", "hal_retz", "mal_hal", "mal_hal_retz"), n_cases = c(sum(mal_hal_rel$all_mal), sum(mal_hal_rel$all_hal),sum(mal_hal_rel$mal_retz), sum(mal_hal_rel$hal_retz), sum(mal_hal_rel$mal_hal), sum(mal_hal_rel$mal_hal_retz)))

# add percent (out of pre + post lines)
line_cat_table <- line_cat_table %>% 
  mutate(perecnt_cases = round((n_cases/length(mal_hal_rel$line_code))*100,1))

## dive deeper into how many people rated what 
# cases where all 3 ratings were present 
mal_hal_retz <- mal_hal_rel %>% filter(mal_hal_retz == 1) %>% select(slide_id, line_code, rating_cols, count_mal, count_hal, count_retzius)
apply(mal_hal_retz[,c("count_mal", "count_hal", "count_retzius")],2,table)

prop.table(table(mal_hal_retz$count_hal, useNA = "always"))
prop.table(table(mal_hal_retz$count_mal, useNA = "always"))
prop.table(table(mal_hal_retz$count_retzius, useNA = "always"))

mal_hal_retz_one <- mal_hal_retz %>% filter(count_hal == 1 | count_retzius == 1 | count_mal == 1) 
length(which(mal_hal_retz_one$count_mal == 1 & mal_hal_retz_one$count_hal > 1 & mal_hal_retz_one$count_retzius > 1)) 
length(which(mal_hal_retz_one$count_hal == 1 & mal_hal_retz_one$count_mal > 1 & mal_hal_retz_one$count_retzius > 1))
length(which(mal_hal_retz_one$count_retzius == 1 & mal_hal_retz_one$count_hal > 1 & mal_hal_retz_one$count_mal > 1)) 
```

### HAL vs. MAL + Retzius Lines

```{r hal vs. mal+retz results, include=FALSE}
## need to collapse ratings into just HAL or other 
hal_other_rel <- mal_hal_rel # create copy (still no dentine bc all are HAL in dentine)

# data pre-check
apply(hal_other_rel[rating_cols],2,table)

# collapse 
hal_other_rel[rating_cols][hal_other_rel[rating_cols] == "mal"] <- "retzius"

# data post-check 
apply(hal_other_rel[rating_cols],2,table)

# lights kappa due to more than 2 raters
hal_other_rel_result <- kappam.light(hal_other_rel[rating_cols])
```

## General Grouping Reliability (Senior/Mid/Junior)  

```{r general category reliability per group, include=FALSE}
## senior ## 
apply(overall_rel[senior],2,table) 

sr_rel_results <- kappam.light(overall_rel[senior])

sr_rel_gwet <- gwet.ac1.raw(overall_rel[senior])

sr_rel_weight <- lights_weighted(raters = senior, weighting = c(0,0.5,1), ratings = mal_hal_rel)

sr_rel_weight_gwet <- gwet.ac1.raw(mal_hal_rel[senior], weights = "linear")

sr_rel_unw_gwet <- gwet.ac1.raw(mal_hal_rel[senior], weights = "unweighted")

## mid ##
apply(overall_rel[mid],2,table)

mid_rel_results <- kappam.light(overall_rel[mid])

mid_rel_gwet <- gwet.ac1.raw(overall_rel[mid])

mid_rel_weight <- lights_weighted(raters = mid, weighting = c(0,0.5,1), ratings = mal_hal_rel)

mid_rel_weight_gwet <- gwet.ac1.raw(mal_hal_rel[mid], weights = "linear")

mid_rel_unw_gwet <- gwet.ac1.raw(mal_hal_rel[mid], weights = "unweighted")

## jr ## 
apply(overall_rel[jr],2,table) # this group consistently has lower yes values

jr_rel_results <- kappam.light(overall_rel[jr])

jr_rel_gwet <- gwet.ac1.raw(overall_rel[jr]) 

jr_rel_weight <- lights_weighted(raters = jr, weighting = c(0,0.5,1), ratings = mal_hal_rel)

jr_rel_weight_gwet <- gwet.ac1.raw(mal_hal_rel[jr], weights = "linear")

jr_rel_unw_gwet <- gwet.ac1.raw(mal_hal_rel[jr], weights = "unweighted")
```

### Seniority Pairings Reliability 

```{r for loop example, include=TRUE, echo=TRUE, eval=TRUE}
sr_mid_table <- data.frame(grouping = character(0), kappa = numeric(0)) # shell dataframe 
for(s in senior){
  for(m in mid){
  x <- kappa2(overall_rel[c(s,m)])
  value <- x$value
  temp <- data.frame(grouping = paste0(s,"+",m), kappa = value)
  sr_mid_table <- rbind(sr_mid_table, temp)
  }}
```

```{r average kappa between general categories, include=FALSE}
## find the kappa between all possible pairing of individual raters between groups and take the average of those 
# calculate each pairings kappa between groups 

## senior x mid ##  
sr_mid_table <- data.frame(grouping = character(0), kappa = numeric(0), gwets = numeric(0), weight_gwets = numeric(0), light_w = numeric(0)) # shell dataframe 
for(s in senior){
  for(m in mid){
  x <- kappa2(overall_rel[c(s,m)])
  kvalue <- x$value
  g <- gwet.ac1.raw(overall_rel[c(s,m)])
  gvalue <- as.numeric(g$est[4])
  gw <- gwet.ac1.raw(mal_hal_rel[,c(s,m)], weights = "linear")
  wg_value <- as.numeric(gw$est[4])
  wl <- lights_weighted(raters = c(s,m), weighting = c(0,0.5,1), ratings = mal_hal_rel)
  temp <- data.frame(grouping = paste0(s,"+",m), kappa = kvalue, gwets = gvalue, weight_gwets = wg_value, light_w = wl)
  sr_mid_table <- rbind(sr_mid_table, temp)
  }}

# AVERAGE of all sr/jr pairings 
sr_mid_avg <- mean(sr_mid_table$kappa) 
sr_mid_avg_g <- mean(sr_mid_table$gwets) 
sr_mid_avg_w <- mean(sr_mid_table$light_w) 
sr_mid_avg_wg <- mean(sr_mid_table$weight_gwets) 

## senior x jr ##
sr_jr_table <- data.frame(grouping = character(0), kappa = numeric(0),gwets = numeric(0), weight_gwets = numeric(0), light_w = numeric(0)) # shell dataframe 

for(s in senior){
  for(j in jr){
  x <- kappa2(overall_rel[c(s,j)])
  kvalue <- x$value
  g <- gwet.ac1.raw(overall_rel[c(s,j)])
  gvalue <- as.numeric(g$est[4])
  gw <- gwet.ac1.raw(mal_hal_rel[,c(s,j)], weights = "linear")
  wg_value <- as.numeric(gw$est[4])
  wl <- lights_weighted(raters = c(s,j), weighting = c(0,0.5,1), ratings = mal_hal_rel)
  temp <- data.frame(grouping = paste0(s,"+",j), kappa = kvalue, gwets = gvalue, weight_gwets = wg_value, light_w = wl)
  sr_jr_table <- rbind(sr_jr_table, temp)
  }}

# AVERAGE of all sr/jr pairings 
sr_jr_avg <- mean(sr_jr_table$kappa) 
sr_jr_avg_g <- mean(sr_jr_table$gwets) 
sr_jr_avg_w <- mean(sr_jr_table$light_w) 
sr_jr_avg_wg <- mean(sr_jr_table$weight_gwets) 


## mid x jr ## 
mid_jr_table <- data.frame(grouping = character(0), kappa = numeric(0), gwets = numeric(0), weight_gwets = numeric(0), light_w = numeric(0)) # shell dataframe 
for(m in mid){
  for(j in jr){
  x <- kappa2(overall_rel[c(m,j)])
  kvalue <- x$value
  g <- gwet.ac1.raw(overall_rel[c(m,j)])
  gvalue <- as.numeric(g$est[4])
  gw <- gwet.ac1.raw(mal_hal_rel[,c(m,j)], weights = "linear")
  wg_value <- as.numeric(gw$est[4])
  wl <- lights_weighted(raters = c(m,j), weighting = c(0,0.5,1), ratings = mal_hal_rel)
  temp <- data.frame(grouping = paste0(m,"+",j), kappa = kvalue, gwets = gvalue, weight_gwets = wg_value, light_w = wl)
  mid_jr_table <- rbind(mid_jr_table, temp)
  }}

# AVERAGE of all mid/jr pairings 
mid_jr_avg <- mean(mid_jr_table$kappa) 
mid_jr_avg_g <- mean(mid_jr_table$gwets)
mid_jr_avg_w <- mean(mid_jr_table$light_w) 
mid_jr_avg_wg <- mean(mid_jr_table$weight_gwets) 
```

## Experience with Deciduous Teeth Groups Reliability  

```{r decid. teeth experience reliability per group, include = FALSE}
## no decid. teeth experience group ##
apply(overall_rel[no_decid],2,table) 

no_decid_rel_results <- kappam.light(overall_rel[no_decid])

no_decid_gwet <- gwet.ac1.raw(overall_rel[no_decid])

no_decid_w <- lights_weighted(raters = no_decid, weighting = c(0,.5,1), ratings = mal_hal_rel)
  
no_decid_wg <- gwet.ac1.raw(mal_hal_rel[no_decid], weights = "linear")

## yes decid. teeth experience group ##

apply(overall_rel[yes_decid],2,table)  

yes_decid_rel_results <- kappam.light(overall_rel[yes_decid])

yes_decid_gwet <- gwet.ac1.raw(overall_rel[yes_decid])

yes_decid_w <- lights_weighted(raters = yes_decid, weighting = c(0,.5,1), ratings = mal_hal_rel)
  
yes_decid_wg <- gwet.ac1.raw(mal_hal_rel[yes_decid], weights = "linear")
```

## Experience with AL Measurement Groups Reliability  

```{r al measure experience reliability per group, include=FALSE}
## no AL measurement experience group ##
apply(overall_rel[no_al],2,table) 

no_al_rel_results <- kappam.light(overall_rel[no_al])

no_al_gwet <- gwet.ac1.raw(overall_rel[no_al])

no_al_w <- lights_weighted(raters = no_al, weighting = c(0,.5,1), ratings = mal_hal_rel)
  
no_al_wg <- gwet.ac1.raw(mal_hal_rel[no_al], weights = "linear")

## yes AL measurement experience group ##
apply(overall_rel[yes_al],2,table) 

yes_al_rel_results <- kappam.light(overall_rel[yes_al])

yes_al_gwet <- gwet.ac1.raw(overall_rel[yes_al])

yes_al_w <- lights_weighted(raters = yes_al, weighting = c(0,.5,1), ratings = mal_hal_rel)
  
yes_al_wg <- gwet.ac1.raw(mal_hal_rel[yes_al], weights = "linear")
```

## Number of Slides Worked on Groups Reliability 
 
```{r number of slides reliability per group, include = FALSE}
## under 200 slides ## 
apply(overall_rel[under_200],2,table) 

under_200_rel_results <- kappam.light(overall_rel[under_200])

under_200_gwet <- gwet.ac1.raw(overall_rel[under_200])

under_200_w <- lights_weighted(raters = under_200, weighting = c(0,0.5,1), ratings = mal_hal_rel)
  
under_200_wg <- gwet.ac1.raw(mal_hal_rel[under_200], weights = "linear")

## 200 -500 slides ##
apply(overall_rel[up_500],2,table) 

up_500_rel_results <- kappam.light(overall_rel[up_500])

up_500_gwet <- gwet.ac1.raw(overall_rel[up_500])

up_500_w <- lights_weighted(raters = up_500, weighting = c(0,0.5,1), ratings = mal_hal_rel)
 
up_500_wg <- gwet.ac1.raw(mal_hal_rel[up_500], weights = "linear")

## 500+ slides ## 
apply(overall_rel[plus_500],2,table)

plus_500_rel_results <- kappam.light(overall_rel[plus_500])

plus_500_gwet <- gwet.ac1.raw(overall_rel[plus_500])

plus_500_w <- lights_weighted(raters = plus_500, weighting = c(0,0.5,1), ratings = mal_hal_rel)
 
plus_500_wg <- gwet.ac1.raw(mal_hal_rel[plus_500], weights = "linear")
```

### Number of Slides Pairing Reliability 
 
```{r average kappa between slide groups, include=FALSE}
# calculate each possible pairing between groups kappa and then take average of all 

## under 200 x 200-500 ##
under_up_table <- data.frame(grouping = character(0), kappa = numeric(0), gwets= numeric(0), weight_gwets = numeric(0), light_w = numeric(0)) # shell dataframe 
for(un in under_200){
  for(up in up_500){
  x <- kappa2(overall_rel[c(un,up)])
  kvalue <- x$value
  g <- gwet.ac1.raw(overall_rel[c(un,up)])
  gvalue <- as.numeric(g$est[4])
  gw <- gwet.ac1.raw(mal_hal_rel[,c(un,up)], weights = "linear")
  wg_value <- as.numeric(gw$est[4])
  wl <- lights_weighted(raters = c(un,up), weighting = c(0,0.5,1), ratings = mal_hal_rel)
  temp <- data.frame(grouping = paste0(un,"+",up), kappa = kvalue, gwets = gvalue, weight_gwets = wg_value, light_w = wl)
  under_up_table <- rbind(under_up_table, temp)
  }}

# AVERAGE of all <200 and 200-500 pairings
under_up_avg <- mean(under_up_table$kappa) 
under_up_avg_g <- mean(under_up_table$gwets) 
under_up_avg_w <- mean(under_up_table$light_w) 
under_up_avg_wg <- mean(under_up_table$weight_gwets) 

## under 200 x 500+ ## 
under_plus_table <- data.frame(grouping = character(0), kappa = numeric(0), gwets = numeric(0), weight_gwets = numeric(0), light_w = numeric(0)) # shell dataframe 
for(un in under_200){
  for(p in plus_500){
  x <- kappa2(overall_rel[c(un,p)])
  kvalue <- x$value
  g <- gwet.ac1.raw(overall_rel[c(un,p)])
  gvalue <- as.numeric(g$est[4])
  gw <- gwet.ac1.raw(mal_hal_rel[,c(un,p)], weights = "linear")
  wg_value <- as.numeric(gw$est[4])
  wl <- lights_weighted(raters = c(un,p), weighting = c(0,0.5,1), ratings = mal_hal_rel)
  temp <- data.frame(grouping = paste0(un,"+",p), kappa = kvalue, gwets = gvalue, weight_gwets = wg_value, light_w = wl)
  under_plus_table <- rbind(under_plus_table, temp)
  }}

# AVERAGE of all <200 and 200-500 pairings
under_plus_avg <- mean(under_plus_table$kappa) 
under_plus_avg_g <- mean(under_plus_table$gwets) 
under_plus_avg_w <- mean(under_plus_table$light_w) 
under_plus_avg_wg <- mean(under_plus_table$weight_gwets) 

## 200-500 x 500+ ## 
up_plus_table <- data.frame(grouping = character(0), kappa = numeric(0), gwets = numeric(0), weight_gwets = numeric(0), light_w = numeric(0)) # shell dataframe 
for(up in up_500){
  for(p in plus_500){
  x <- kappa2(overall_rel[c(up,p)])
  kvalue <- x$value
  g <- gwet.ac1.raw(overall_rel[c(up,p)])
  gvalue <- as.numeric(g$est[4])
  gw <- gwet.ac1.raw(mal_hal_rel[,c(up,p)], weights = "linear")
  wg_value <- as.numeric(gw$est[4])
  wl <- lights_weighted(raters = c(un,p), weighting = c(0,0.5,1), ratings = mal_hal_rel)
  temp <- data.frame(grouping = paste0(up,"+",p), kappa = kvalue, gwets = gvalue, weight_gwets = wg_value, light_w = wl)
  up_plus_table <- rbind(up_plus_table, temp)
  }}

# AVERAGE of all <200 and 200-500 pairings
up_plus_avg <- mean(up_plus_table$kappa) 
up_plus_avg_g <- mean(up_plus_table$gwets) 
up_plus_avg_w <- mean(up_plus_table$light_w) 
up_plus_avg_wg <- mean(up_plus_table$weight_gwets) 
```

```{r check total markings per rater, include=FALSE}
rater_counts <- data.frame(rater = c(senior, mid, jr),
                seniority = c("senior","senior","mid","mid","mid","junior","junior","junior"),
                al_experience = c("no","yes","yes","yes","yes","no","no","no"),
                deciduous_teeth_experience = c("no","yes","yes","yes","no","no","no","no"),
                total_slides = c("500+","500+","200-500","200-500","<200","<200","<200","200-500"))

# shell
shell <- data.frame(rater = NA, n_retzius = NA, n_mal = NA, n_hal = NA)
for (i in rating_cols){
  retz <- length(which(mal_hal_rel[,i] == "retzius"))
  mal <- length(which(mal_hal_rel[,i] == "mal"))
  hal <- length(which(mal_hal_rel[,i] == "hal"))
  temp <- data.frame(rater = i, n_retzius = retz, n_mal = mal, n_hal = hal)
  shell <- rbind(shell, temp)
}

shell <- shell %>% filter(!is.na(rater))
rater_counts <- left_join(rater_counts, shell, by = "rater")
```


